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<^ ABSTRACT 

Radiative transfer (RT) problems in which the source function includes a 
scattering-like integral are typical two-points boundary problems. Their solution 
1. via differential equations implies to make hypotheses on the solution itself, namely 

j_| ■ the specific intensity /(r; n) of the radiation field. On the contrary, integral 

^ methods require to make hypotheses on the source function S{t). It looks of 

course more reasonable to make hypotheses on the latter because one can expect 
that the run of S{t) with depth be smoother than that of /(r; n). 

In previous works we assumed a piece-wise parabolic approximation for the 
source function, which warrants the continuity of S{t) and its first derivative at 
each depth point. Here we impose the continuity of the second derivative S"{t). 
In other words, we adopt a cubic spline representation to the source function, 
which highly stabilize the numerical processes. 
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1. Introduction 

Some years ago we proposed a new algorithm, the Imphcit Integral Method (IIM), to 
solving those radiative transfer problems in which the specific source functions (one for each 
frequency and direction pair) depend linearly on the radiation field via a single quantity 
independent of both frequency and direction. In the paradigm instance of radiative transfer 
through an ideal medium formed by atoms with only two energy levels (Two-Level Atom 
model), this quantity is the integral over frequencies of the mean specific intensity of the 
radiation field, weighted with the spectral profile. (See Simonneau and Crivellari, 1993, 
hereinafter Paper I.) 

Because it is independent of both frequency and direction, such a quantity constitutes 
a single scalar coupling for all the specific RT equations, and can be chosen in a natural way 

as the protagonist variable for the numerical solution of the RT problem. This choice is the 
distinctive and essential feature of our IIM: to work with a quantity which is independent 
of both frequency and direction brings about that the method does not require to store and 
invert huge matrices like in the customary numerical algorithms employed in RT problems. 
We have already remarked in Paper I that our algorithm is a mere phenomenological repre- 
sentation of the actual physical process. Because of that and due to the lack of a matricial 
structure, the advantages of the IIM in terms of reliability, accuracy and robustness sould 
be self-evident, as well as the conspicous saving of both computational time and memory 
storage it makes possible. 

The aforesaid advantages suggested us the possibility to employ the IIM also in the 
computation of stellar atmospheres models, where we must solve many (some hundreds) RT 
equations, one for each frequency. The source function of each specific RT equation is here 
the weighted mean of a term that includes the mean specific intensity of the radiation field 
through a scattering-like integral with a thermal contribution given by the Planck function 
By{T). The paradigm problem of the self-consistent temperature correction when computing 
stellar atmosphere models was considered in Crivellari and Simonneau (1994). 

First of all, we must recognize that the geometrical structure of the system, that is the 
sequence of the discrete atmospheric layers, must be necessarily the same for all the frequen- 
cies. But we must also recognize that for any given frequency some layers do not contribute 
to the formation of the spectrum. They do not take an effective part in the radiative transfer 
process because either they are exceedingly transparent {i.e. exp{— Ar,^} ~ 1) or they corre- 
spond to optically very deep regions (i.e. exp{— Ar,^} ~ 0). The layers intermediate between 
the above two groups constitute the specific spectral formation region. However all the layers 
of the structure must be taken into account in the numerical algorithm, irrespectively of the 
frequency considered. Yet, due to the dramatic difference among the values of the opacity 
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with frequency, different spectral intervals form in very different geometrical regions. That 
compels us to divide the atmosphere into very many layers in order to cover properly all 
the spectral formation intervals. On the other hand, it is matter of the run with depth of 
the data that are common to radiative transfer at all the frequencies. As an exemple, given 
a temperature distribution on the discrete atmospheric layers, the variation with depth of 
the numerical values of the Planck function By{T), i.e. the monochromatic thermal sources, 
may vary enormously frequency by frequency. For instance, in the case of a solar-like star, 
from the bottom to the top of the atmosphere Bi,{T) varies by a factor of the order of 10^ 
for frequencies in the visible part of the spectrum, while this factor can be of the order of 
10^^ for frequencies in the range of Lyman a. Therefore a set of depth points suitable for a 
good description of the mathematical behaviour of the source function at some frequencies 
cannot be adequate at other frequencies. Again very many common discrete depth points 
are necessary in order to provide a proper distribution of the data for the adequate treatment 
of each monochromatic RT equation. 

The foregoing requirements make it impossible in the practice to replace derivatives by 
finite differences, as in the outermost layers the optical thickness is almost zero for many 
frequencies. The use of integral methods may seem to be the only advisable way out, but the 
very large number of discrete optical depth points, necessary to warrant the proper treatment 
of the RT process at all the frequencies, does not advise to employ global integral methods, 
too. 

We can get rid of the difficulties brought about by the introduction of very many layers 
on the one hand by employing our IIM, which allows us to take into consideration as many 
geometrical depths as necessary because, as already said, it does not require the storage and 
inversion of huge matrices. Moreover, on the other hand, we can introduce a better mathe- 
matical representation of each monochromatic source function Si,{t^) in order to account for 
the possible rapid variation of both the branching parameter e^, (see eq. [2] later) and B^{T) 
with respect to each specific optical depth Tj,. In such a way we can optimize the treatment 
of all the individual frequencies. 

In the original formulation of the IIM (see the above references) we considered models 
that comprised 150-200 discrete layers between the surface and the bottom of the atmosphere. 
Inside each of them we approximated each specific source function Su{t^) by an arc of 
parabola and imposed the continuity of S,y{T„) and Sl{T,y) at all the NL dividing points. This 
piece-wise parabolic approximation yielded excellent results in many cases (see the above 
references). However under extreme conditions, for instance in the case of sudden variations 
of the thermal sources {e.g. at Lyman a frequencies in cool stars), such an approximation 
can introduce numerical instabilities that spoil the computation of the model. 
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To impose also the continuity of the second derivative of Si,{t^) at all the NL dividing 
points can remove the foregoing instabilities. Consequently we propose here a cubic spline 
model for each specific source function. In some way this model constitutes a regularization 
of the process to computing the values of the source functions. The formalism of the cubic 
spline approximation (namely a two-point boundary value problem developed to interpolate 
among the NL explicitly known values of a given function) can be employed in the present 
case although the NL values of SyiT^) are yet unknown. 

To employ the cubic spline approach in order to describe the behaviour of the source 
function in typical RT problems, where a scattering term appears in the source function, 
is the best (may be the unique) correct choice for both theoretical and numerical reasons. 
A theoretical reason is brought about by the non-local nature of the problem: the specific 
intensities and consequently the source function at a given depth point depend via the RT 
process on the values of the source function at all the other points of the system. Thus the 
numerical values of the source function must be computed simultaneously at all the depth 
points. Therefore such a non-local character of the physical problem must be represented by 
means of a non-local mathematical structure. Also the derivatives of the source function at 
any depth point must be formulated as a linear relation including the implicit values of the 
source function at all the depth points, not only as a linear relation of the implicit values of 
the source function at each triad of consecutive depth points. 

The pratical reason is for the sake of the stability of the computational algorithm. 
The cubic spline model minimizes the strain energy integral, that is the integral of the 
squared values of the values of the second derivative of the protagonist function, namely 
of the variation of its curvature - i.e. the oscillations. (See, e.g., Rivlin, 1981.) That 
is, the use of the cubic spline approximation to the source function minimizes the risk of 
destabilizing oscillations. The cubic spline representation constitutes by itself another two- 
points boundary problem: we must know the values of S{t) at all the depth points in order 
to compute S'ij) and S"{t) at each depth point. Nevertheless this apparent drawback 
turns out to be on the contrary an advantage, because we can carry on simultaneously both 
two-points boundary problems: the RT problem and the cubic spline interpolation. 

Prom the algorithmical stand point, the kernel of the original IIM is a forward- elimination 
scheme that links the so far unknown values of the source function at each pair of consecu- 
tive optical depth points (r^, Tl+i) by mean of a linear relation with known coefficients. The 
latter are determined by taking into account the RT equations that describe layer by layer 
the propagation of both the downgoing and the upgoing specific intensities. Now we realized 
that, by using the cubic spline formalism the same forward-elimination scheme can also be 
employed to link the unknown values of the second derivatives of the source functions, again 
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by means of a linear relation. 

Once attained the deepest optical depth point ttvl at the end of the forward-elimination, 
we can impose the bottom boundary condition (eq. [5] later on) to both the RT process and 
the cubic spline chain; in other words we can close the linear relation between Si,{tnl-i) 
and S,^{tnl) on the one hand, between S"{tnl-i) and S'1{tnl) on the other. This al- 
lows us to recover the numerical values of the source functions and their second deriva- 
tives at the bottom, as well as those of the set of the incident upgoing specific intensities 
{I^{Ti\fL, fij), J = 1,ND}. Then, in a succesive back-substitution scheme, we are in a posi- 
tion to compute at each depth point the numerical values of the source functions and their 
second derivatives by using the above linear relations, whose coefficients have been stored 
during the previous forward-elimination. 

Thanks to that we have at hand a unique algorithm to solve each specific RT problem 
under the imposed constraint that the specific source functions as well as their first and 
second derivatives be continuous at all the NL points of the grid chosen for the geometrical 
representation of the stellar atmosphere. In such a way we can get rid of the instabilities 
that may arise in the case of extreme variations of the source functions without paying any 
extra computational cost. 



2. The mathematical background 

For the sake of an easier presentation of the new more precise version of the IIM an- 
nounced in Section 1, we will consider the simplest instance that yet contains all the dif- 
ficulties intrinsic to RT astrophysical problems, namely the transport of monochromatic 
radiation through a plane-parallel medium in which matter particles can scatter, absorb and 
emit photons. In the previous works above quoted the original formulation was applied to 
much more general instances. The version presented here can be easily applied to such cases. 

Following the customary notation, the RT equations that describe the evolution of the 
upgoing intensities /"^(r, //) and the downgoing intensities /~(t, //) are 

±^,±I^{r,^)= I^{t,^) - S{t) , (1) 

where r denotes the optical depth and fi is the cosine of the angle formed by the direction 
of propagation with the perpendicular to the plane-parallel layers (// = cos 0,0 < /i < 1) 

The source function is a weighted mean between the thermal source B{t) and the 
mean intensity J(t), namely 
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S{t) = eB{T) + (l-£)J(r). (2) 

The branching parameter e = 6{t) is the ratio of the absorption coefficient to the 
total opacity {i.e. the sum of the absorption and the scattering coefficient). The latter 
defines the scale of the optical depth r; (1 — e) is customarily called the albedo. In terms of 
the upgoing and the downgoing intensities the mean intensity is given by 

Ar) = [\l^{r,i,) + /-(T,/x)]d/.. (3) 

^0 

The integral in eq. (3) is representative of any scattering integral, which may be different 
for the application of the IIM to different instances. 

In the discrete ordinates approximation the integral in eq. (3) is replaced by the sum 
of the intensities corresponding to a finite number of ND directions. Then 

ND 

J{t) «^ J]wj[7+(t,/xj) + /-(t,//j)] . (4) 
j=i 

For most RT problems in plane-parallel geometry (at least for stellar atmosphere 
models computations) a five-points Gauss division of the interval < // < 1 is more than 
enough. The Wj's are the corresponding integration weights. 

The numerical solution requires the discretization of the optical depth variable r, 
too. The stellar atmosphere must be sliced into a set of NL plane-parallel horizontal layers, 
divided by the set of NL + 1 optical depths points {TQ,ri,T2...rNi}. The value Tq = 
corresponds to the surface and tj^l to the bottom of the atmosphere. The computation of a 
fairly good model require that NL be of the order of two hundred. 

The values of the incident intensities onto the top surface, i.e. the downgoing intensi- 
ties I~{tq, fij), and those of the incident intensities onto the bottom surface, i.e. the upgoing 
intensities I~^{tnl, fJ^j), must be known; they are data of the RT problem. In the case of 
a stellar atmosphere /~(ro,/i) is usually assumed to be zero, that is there is not radiation 
incident onto the stellar surface. We will show later that the method can equally work also 
under more general conditions. For the upgoing intensities at the bottom of the atmosphere 
we can assume that the diffusion approximation holds valid, that is 



(5) 
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which is brought about by the cubic polynomial behaviour of S{t) at depths immediately 
greater than Tp^^. These two families of boundary conditions are sufficient to ensure that 
the RT problem is self-consistent. 

The link between the values of the specific intensities at any pair of consecutive 
optical depth points (tl, tl+i), namely any single link of the whole RT chain, is given by the 
corresponding RT equations in the integral form, that is 



/+(tl,/^j) = /+(Ti+i,//^)exp(-^) + Sit)eM-^-^)dt (6) 

and 

I~{tl+i,I1j) = /"(TL,/xj)exp(-^^) + / S{t)cxp{-^^ — -) dt , (7) 

where Atl = tl+i — tl- Equations (6) and (7) are the straightforward representation of the 
RT process. 

At the surface (i.e. for Tq = 0) the set of values {I^{0, fij), J = 1,ND} are the 
initial conditions for the inward RT problem, while the set {/+(0,/ij), J = 1,ND} is the 
solution of the outward RT problem, i.e. the emergent intensities. At the bottom the 
set {I~^{tnl, fJ-j): J = 1,-^-0} yields the upgoing initial conditions (cf. eq. [5]); the set 
{I~{tnl: l^j): J — N D} IS the result of the inward RT process. 

Let us now turn our attention on the cubic spline approximation to the source function 
S{t). That is, we will assume a cubic polynomial approximation inside each particular 
interval {tl,tl+i), defined by two consecutive optical depth points, as the single link of the 
spline chain. Anyone of these arcs of cubic is uniquely determined by the values of the source 
function and those of its second derivative at the end points (knots) of each interval. 

To impose the continuity of the source function as well as that of its first and second 
derivative at the end points of each interval {tl,tl+i) leads to the cubic spline condition 



^ S"(n_,) + + ^) S'VJ + ^ S"(r.«) , (8) 
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where Atl = tl — tl-i and Atl+i = tl+i — tl- Likewise, as a consequence of the cubic 
behaviour of S{t) between tl and r^+i, it will hold that 

S'{tl) = ^ [^(TL+i) - S{t,)] - ^S"{t,) - ^S"{t,^,) , (9) 



(10) 



and 



(11) 



By means of eq. (8) we are in a position to join the neighbouring links of the spline chain 
while ensuring the required continuity at the knots. 

Like for the RT chain, also for the cubic spline chain wc need two bundary conditions. 
Customarily these are S"{to) = S"{ti) and S'^tj^l) = S"{tj^l_i). In the present study we 
assume that S"{to) = S"{ti) at the surface, namely that the first arc of the spline chain is a 
parabola. On the contrary, the boundary condition for the spline chain at the bottom must be 
consistent with the diffusion approximation for the incident upgoing intensities /"^(tjvl, /ij), 
given by eq. (5), which is a consequence of having assumed also a cubic polynomial behaviour 
for S{t) at depths greater than tnl- This condition is in agreement with the cubic polynomial 
behaviour of S{t) inside the last layer (tjvl-Ij ^Afi)- Hence we cannot introduce now a 
different approach to S{t). However we can derive the formal value of the first derivative 
S'{t) from eq. (2), that is 

S'{r) = [e{r)B{r)] ' + [1 - e{r)] ' J{r) + [1 - e{r)] J'{r) , (12) 

where 

Equation (12), evaluated at the deepest optical depth point tnl will then yield the required 
lower boundary condition, as will be shown later. 

Let us get back now to eq.s (6) and (7). For any interval {tl,tl+i) the arc of cubic 
approximating to S{t) is given by 
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S{t) = S{tl) + S'{tl){t - tl) + \s"{tl){t - r^f + \s"\tl){t - tl)\ (14) 

By replacing eq. (14) in eq.s (6) and (7), and taking into account eq.s (9) through (11), we 
get eventually 



I^{TL,idj) = /+(rL+i,/xj)exp-(^^) + 
WSi( J) S{tl) + WS2( J) S{tl+i) + wdi( J) S"{tl) + Wd2( J) S"{tl+i) (15) 

and 



I (tl+i,//j) = I {TL,f^j)e^-{ )+ 

WS2( J) S{tl) + WSi(J) S{tl+i) + wd2(J) S"{tl) + wdi(J) S"{tl+i). (16) 

The quadrature weights wsx(t/), ws2(t/), wd]^(t/) and wd2(t/) are computed straightforwardly 
by taking into account eq.s (9) through (11) to yield 

wsi(J) = [l-l] + ie-^ (17) 

WS2(J) = ^ - [1 + ^] e"'' (18) 

wdi(J) = [(1 - ^ - ^) - (^ - ^) e-^] , (19) 

wd2(J) = [(-^ + 1) - (1 + ^ + 1) e-'] , (20) 

where S — Atl / ixj. 

Sometimes, when 5 « 1, for sake of numerical percision it may be necessary to recast 
the foregoing weights into the form 
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wsi(J) = -6 - -6' + -S' - —5^ + —S' + ^—5^ , (21) 

^ ^ 2 6 24 120 720 5040 40320 ' ^ ^ 

wd,(.) = - (1.= - + ±^S^ - Jj.' + . (23) 

wd.(.) = - (1.3 - + - + 3!^.') . (24) 

To conclude, eq.s (15) and (16) together with eq.s (17) through (20) allow us to write 

explicitly for each direction /ij the relations between I^{TL,iJij) and /'''(Ti+i, /xj) on the one 
hand, between I~{tl^i, fij) and I~{tl,^j) on the other. These relations are linear functions 
of the unknown values of S{tl), S{tl^i), S"{tl) and S"{tl^i), which will play a protagonist 
role in the numerical algorithm. The cubic spline condition, given by eq. (8), impose a 
further relation between S{tl) and S"{tl) at each knot tl- 

We recall that for any frequency the specific source function is approximated by 
an arc of cubic inside each interval {tljTl^i). Therefore in the layers deeper than the 
corresponding spectral formation region, where exp{— Ar} is pratically null, the form of the 
weights wsi(J), ws2(J), wdi(J) and wd2(J) given by eq.s (17) through (20) warrants that 
the intensities I~^{t,iij) recover there the form of eq. (5), originally assigned at the bottom 
of the atmosphere {i.e. at tjvl). That is to say, the boundary condition, initially assigned at 
the bottom, is transported up to the end of the spectral formation region, keeping its form in 
a natiual way. On the other hand, in the outer layers beyond the region of formation,where 
exp{— Ar} approaches unity, eq.s (6) and (7) warrant that /"""(r, /xj) and /^(r, /xj) keep 
constant. Thus, albeit the total number NL of layers exceed that required by the proper 
physical treatment of the formation region for each single frequency, such an excess docs 
not affect the numerical computation of the protagonist variables. That is to say, frequency 
by frequency the effective transport of the specific intensities is performed in a natural way 
inside its own region of formation, provided that care has be taken to select the geometrical 
width of the stellar atmosphere system so that, as already stressed in the Introduction, the 
former include the region of formation for all the frequencies. 
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We have then at hand all the mathematical tools that will allow us to solve the global 
RT problem in the same way as in the original IIM scheme (see Paper I) . 

However only to warrant the continuity of the two first derivatives of the source func- 
tion is not enough to avoid the occurence of instabilities. As in the cubic spline fundamental 
equation (eq. [8]) the protagonist variables are the function itself and its second deriva- 
tive (both tied through their values at any set of three consecutive points), also in the RT 
elimination scheme the source function and its second derivative must be the protagonist 
variables. 

In a previous attempt we formulated the equations (15) and (16), which describe the 
propagation of the upgoing and downgoing intensities, in terms of S{rL), ^(r/^+i), S'{rL) and 
S'{tl+i) after the elimination of S"{tl) and S"{tl^i) given as functions of S{tl), Siji^i), 
S'{tl) and S'{tl^i) thanks to the cubic behaviour of S{t). In the actual version we describe 
the propagation of the aforesaid intensities by means of S{tl), S{tl+i), S"{tl) and S"{tl+i)i 
again by means of the cubic behaviour of S{t). From the mathematical point of view both 
representations should yield the same results, but from the numerical standpoint it looks 
much better to work directly with the second derivatives S"{tl) and S"{tl^-i), because the 
fundamental equation (8), that links the sequence of succesive layers in the cubic spline 
scheme, requires the variables S{t) and S"{t). 

In the present formulation of the propagation equations (15) and (16) the integration 
weights wsi(J) and WS2( J), given by eq.s (17) and (18), account strictly for the hnear piece- 
wise approximation to any monochromatic source function Sv{t). The remaining weights, 
wdi(J) and wd2(J), account for the deviation from the linear behaviour, either parabolic 
or cubic. Whenever wdi(J) and wd2(J) take on small values, the linear approximation is 
more than enough. This is the case in the outermost layers, where it holds that 5 < 1 and 
5^ « 1; the linear approximation is automatically recovered, as only the weights wsi(j) 
and WS2(J) account for the variation of the source function in optically thin layers. That 
is to say, in the practice only S{ti^) and SXr^+i) take part in the elimination scheme. In 
other words, the effects of a non-linear behaviour play the role of a perturbation of the linear 
behaviour. 

In the original formulation of the IIM (Paper I), we employed the variables S{tl), 
S{tl+i), S'{tl) and S"(rL+i), together with the corresponding integration weights, in order 
to describe the propagation of the upgoing and downgoing intensities between any pair of 
optical depth points tl and r^+i. Whatever their behaviour (linear, quadratic or cubic), all 
the four variables and the relevant integration weights took an active part in the elimination 
scheme, both from the theoretical and the numerical standpoint. This can have been at the 
origin of the instabilities that showed up, above all in the regions of small optical depth. 
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The actual version of the IIM, due to the above mentioned reasons, results certainly more 
rehable. 

3. The Forward-Elimination/Back-Substitution scheme 

As already said, we will work with a set of fundamental variables whose values are un- 
known: the upgoing and downgoing specific intensities /^(tl,//j), the corresponding source 
functions S{tl) and their second derivatives S"{tl)- The major aim of this section is to 
derive linear relations among the values of the foregoing fundamental variables at the two 
consecutive optical depth points tl and r^+i that delimitate each of the layer {jl^tl+i) 
succesively under study. The coefficients of these relations are easily computed, and will be 
denoted in the following by bold face symbols. 

3.1. The algorithmic representation of the upper boundairy conditions 

Wc start necessarily with only one half of the data of the problem, namely the set of the 
downgoing intensities incident onto the upper boundary layer at Tq = 0, i.e. {/"(tq, /Xj), J — 
1, ND}, that we will write in its most general form as 



r{To,fXj) = cmO(J) + cmsl(J) S{to) + cms2(J) S{ti) + 

ND 

cmdsl(J) S"{to) +cmds2(J) S"{ti) + ^ R( J, J') /+(to,/^j') . (25) 

j'=i 

The coefficient cmO(J) accounts for the numerical value of the incident intensity /^(ro,/ij), 
which is usually null. The reflexion matrix R( J, J') takes into account the possible ef- 
fects of backscattering outside the stellar surface. Under usual conditions it holds that 
also R(J, J') = 0. On physical grounds it is hard to justify the dependence of I~{to,ij,j) 
on the values of the source function and its second derivative at points tq and ti through 
the coefficients cmsl( J), cms2( J), cmdsl( J) and cmds2(J). It is rather an algorithmical 
requirement, as these coefficients allow us to link linearly the values of the protagonist vari- 
ables between two consecutive optical depth points. Consistently with the upper boundary 
conditions, the latter coefficients have to be set equal to zero. 

At the end of the treatment of radiative transfer in the first layer (as well as in the 
succesive ones) some of these coefficients will take on values different from zero. These 
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new values can overrun the previous memory storage, because the current relation for the 
downgoing intensities at Tq = will not be necessary any longer. 

Inside the forward-elimination scheme for the RT process we must propagate not only 
the upgoing and downgoing specific intensities (which brings about the propagation of the 
source function as defined by eq.s (2) and (3)), but also the second derivative S"{t) of the 
source function in the cubic spline scheme. 

As already said, we assume that in the first layer (tq, Ti) the source function S{t) can 
be approximated by an arc of parabola, which imphes that S"{to) — S"{ti). This condition 
will be included in the coefficients of the relation 

S"{to) = cdsO + cdsl S{to) + cds2 S{ti) + 

ND 

cddsl 5"(to) + cdds2 5"(Ti) + J] cdi(J) /+(ro,/ij) , (26) 

j=i 

where the values of S{tq), S{ti), S"{tq), S"{ti) and the set {I^{tq, jij), J = 1,ND} are un- 
known. In order to fulfill the above boundary condition, all the coefficients in cq. (26) must 
be equal to zero, excepted cdss2 that must be set equal to one. To express here S"{to) as a 
function of S{to), S{ti), S"{to) itself, S"{ri) and the set {/+(ro,/ij), J = l,ND} is just for 
algorithmical ease. When convenient, we will solve for S{to) - and for S"{to) - in terms of 
5(ri), 5"(ti) and 

3.2. The layer by layer elimination 

We are going to show here how the treatment of the first layer (ro,ri), labelled by L = 1, 
will yield the coefficients of the relation 

S{to) = cbsO(l) + 

ND 

cbss(l)5(Ti) + Cbsd(l) 5"(ti) + J2 Cbsi(l,j) /+(Ti,/ij) (27) 

J=l 

and those of the relation 

S"{to) = cbdO(l) + 

ND 

cbds(l) S{ti) + cbdd(l) S"{ti) + ^ cbdi(l, J) /+(ri,/xj) . (28) 

j=i 
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These coefficients will be stored in order to compute S{to) and S"{to) in the succesive back- 
substitution process, once the values of S{ti), S"{ti) as well as the set {I~^{ti, /ij), J = 
1,ND} have been determined. The above relations link any pair of succesive layers. As 
already said, the determination of these relations constitutes the aim of this section. 

In parallel we are going to show also how to recover the initial conditions for /~(ti, /ij) 
and S"{ti), i.e. the values of the coefficients of the relations equivalent to eq.s (25) and (26), 
now for Ti. 

Let us detail our foregoing purpose. At the beginning of the study of each succesive 
layer - here the first one - we must consider the implicit computation of the corresponding 
source function at the upper limiting optical depth, here Tq. The form of the incident 
downgoing intensities at tq, given by eq. (25), together with the implicit values of the set 
{I~^{tq, fij), J = 1, ND} allow us to compute from eq. (4) the coefficients of a linear relation 
among J(ro) and S{to), S{ti), S"{to), S"{ti) and the set {I'^{to, /ij), J = 1,ND}. Then eq. 
(2), where s{to) and -B(to) are given, will yield the coefficients of the linear relation 

S{to) = csO + cssl S'(to) + css2 S{ri) + 

ND 

csdsl S"{to) + csds2 S"{ti) + ^ csi(J) I^{ro,nj) ■ (29) 

j=i 

where we have not solved for S{to) again for the sake of algorithmical ease. 

We compute now for each direction /ij the quadrature weights wsl(J), ws2(J), 
wdl(J) and wd2(J) according to eq.s (17) through (20) - or alternatively eq.s (21) through 
(24) - for Ati = Ti — Tq. These weights allow us an implicit quadrature of the source function 
in the description of the propagation of the upgoing intensities from /+(ti, /Xj) to I~^{to, /Ij), 
and later of the downgoing intensities from I~{to,ij,j) to I~{ti,ij,j). 

At this point we can introduce the implicit form for I^{tq,ixj) in terms of S'(ro), 
5(ri), 5"(ro), S"{ti) and the set {/+(ri,/ij), J = 1,ND}, given by eq. (15), in eq. (25) for 
I~{tq, fij), which is the initial condition for the study of the layer (tq, Ti). By re- arrangement 
of the coefficients we can write 

/~(to,//j) = cmO(J) + cmsl(J) S(ro) + cms2(J) S(ri) + 

ND 

cmdsl(J) S"{to) +cmds2(J) S"{ti) + ^ R( J, J') I+{ti,iij>) . (30) 

j'=i 

for any direction /ij. These new values of the coefficients can overrun the memory places of 
the pervios ones, corresponding to the intial condition given by eq. (25). 
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We repeat the same exercice, namely to employ eq. (15) inside both the functional 
form for S{to), given by eq. (29), and that for S"{tq), given by eq. (26), in order to recover 
the previous form for both of them, but now as a function of the upgoing intensities at Ti 
insted of Tq, hence with different coeffficients. That is 

S{to) = csO + cssl S{to) + css2 S{ti) + 

ND 

csdsl S"{ro) + csds2 S"{ti) + ^^^i'^) ^^i^uf^j) ■ (31) 

j=i 

and 



S"{to) = cdsO + cdsl S{to) + cds2 S{n) + 

ND 

cddsl S"{ro) + cdds2 ,5"(ti) + ^ cdi(J) /+(Ti,Atj) . (32) 

j=i 

Now, just by solving for S{to) and S"{to) we obtain the coefficients of the relations (27) and 
(28), earlier announced at the beginning of Section 3.2. These coefficients must be stored 
for further use. In such a way we have achieved part of out goal. 

At this point let us describe the propagation of the downgoing intensities from /^(tq, /xj), 
given by eq. (30), to /~(ri,//j) according to eq. (16). By re-arrangement of the coefficients 
we get the new values corresponding to the relation 

I~{n,iJ,j) = cmO(J) + cmsl(J) 5'(ro) + cms2(J) ^(ri) + 

ND 

cmdsl(J) 5"(ro) + cmds2(J) 5"(ri) + ^ R{J, J') I+{n, (33) 

j'=i 

for all the directions ^j. Again these new values of the coefficients can overrun the previous 
ones, corresponding to eq. (30). 

If we introduce the foregoing eq.s (27) and (28), whose coefficients we have just 
computed, in the functional form of /~(ti,;Uj) given by eq. (33), by re-arrangement of the 
previous coefficients we derive the new ones for the relation 

ND 

cmO(J) + cmsl(J) S{ti) + cmdsl(J) S"{ti) +^ R( J, J') /+(ri,/^j/) , (34) 

j'=i 
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which we will cast into the form required by eq. (25) by setting equal to zero the coefficients 
cms2(J) and cmds2(J). We have thus determined the coefficients of the linear relation 
required as the initial condition at ri, that will be necessary to study the propagation of the 
downgoing intensities in the succesive layer (ti,T2). 

We have still to determine the initial condition for the propagation of 5"'(r), that is 
to say a linear relation like cq. (26), now for Ti. It is matter of recovering the functional 
form of S"{ti) in ordcrt to start the study of the spline chain in the layer (ri, T2). We have at 
hand the fundamental relation for the cubic spline, namely eq. (8) that links linearly S"{tq), 
S"{ti) and S"{t2) with S{to), S{ti) and ^(rs). 

By introducing in eq. (8) the formal expressions for <S'(to) and ^"'(to), given by eq.s 
(27) and (28), we get easily the coefficients of the equation 

S"{ti) = cdsO + cdsl S{ti) + cds2 ^(ra) + 

ND 

cddsl5"(ri) + cdds2,S"(T2) + J2 "^^^('^) ^^(^i'/^^) ' (35) 

j=i 

akin to eq. (26), the bootstrap at tq, that was the initial condition to studying the layer 
(to,ti). Equation (35), together with eq. (33) that is the initial condition for the treatment 
of radiative transfer, will allow us to repeat the foregoing procedure for the layer (ri,T2). 
This scheme is then iterated layer by layer till the bottom of the atmosphere. 

3.3. The solution at the bottom and the Back-Substitution 

At the end of the forward-elimination scheme we have at hand the full set of coefficients 
of eq.s (27) and (28) for each optical depth of the set {to,ti, ...,tnl-i}- The expUcit values 
of -S'(Ti) and S"{tl) as well as those of the set of the outgoing intensities {I~^{tl, Hj), J — 
1, ND} have now to be computed in the back-substitution scheme. 

For the sake of a more clear exposition of the mathematical solution at the bottom 
we will rewrite eq.s (27) and (28) for tatl, that is 

S(tnl-i) = cbsO(A^L - 1) + 

ND 

chss{N L - 1) S{rNL) + chsd{N L - 1) S"{rNL) + J] cbsi(iVL - 1, J) /+(r^L,/ij) (36) 

j=i 
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and 

S"{tnl-i) = cbdO(7VL - 1) + 

ND 

cbds(7VL-l) S{tnl) + cbdd(7VL-l) S"{tnl) + cbdi(iVL-l, J) /+(rjvL,/^j) • (37) 

j=i 

Also, at the end of the forward-ehmination scheme, the current values of the coefficients of 
the equation for I~{tnl, l^j), that is 

I~{rNL-,IJ'j) = cmO(J) + cmsl(J) S{tnl) + cms2(J) S{tnl+i) + 

cmdsl(J) S"{tnl) +cmds2(J) S"{tnl+i) + ^ R( J, J') /+(riVL,/xjO , (38) 

J' 

are still stored in the scratch memory. For the sake of a homogeneous algorithm we had kept 
the dependence on S{r]siL+i) and S"{t]^l+i) through the coefficients cms2( J) and cmds2( J). 
But these coefficients are null so that S{t]^l+i) and S"{t]^lj^i) do not play any active role. 
The same algorithmical requirement compelled us to introduce the dummy supplementary 
optical depth ttvl+i- 

At this point we can apply the lower boundary condition for the radiative transfer, 
i.e. the formal expression for I^{tnl, fJ'j) given by cq. (5). If wc replace this expression 
in the previous cq.s (36), (37) and (38), by re- arrangement of terms we obtain the explicit 
values of the coefficients of the two linear relations for S{tni-i) and S"{tnl-i) as a function 
of S{T]\fL), S'{tnl), S"{Tj\fL) and S"'{tnl), that we will write as 

S{tnL-i) = Cn{S{TNL),S'{TNL),S"{Tr,L),S"'{TNL)} (39) 

and 

S"{rNL-l) = Cn{S{Tr,L):S'{Tr,L):S"{Tr,L):S"'{rj,L)} . (40) 

Likewise, if we take into account the aforesaid expression for I~{tnl, IJ'j), whose coefficients 
cms2(J) and cmds2(J) are null, we get also the coefficients of the hnear relations 

I-{TNL,fij) = Cn{S{TNL),S'{TNL),S"{TNL),S"'{TNL)} (41) 

for each direction /ij. 

In the forward-elimination, at the beginning of the study of each layer {tl,tl+i), we 
have formally computed the mean intensity J{tl) and the corresponding source function 
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S{tl) at the upper optical depth tl- That is to say, we have not yet used the relation given 
by eq. (2) at the last optical depth tatl- Let us do it here. 

Equation (5) for I^{tnl, I^-j) and (41) for I~{tnl, Hj) allow us to compute via eq. (4) 
the coefficients of a linear relation like 

J{tnl) = Cn{S{TNL),S\TNL),S'\TNL),S"\TNL)} (42) 

Now thanks to eq. (13) we can compute also the coefficients of the linear relation 

J'{tnl) = Cn{S{TNL),S\TNL),S'\TNL),S"\TNL)} . (43) 

By means of eq.s (2) and (12) we can eventually derive the explicit coefficients of the linear 
relations 

S^Tnl) = Cn{S{TNL),S\TNL),S"{TNL),S"\TNL)} (44) 

and 

S\tnl) = Cn{S{TNL),S\TNL),S"{TNL),S"\TNL)} . (45) 

The two latter relations are the independent conditions to close both the radiative transfer 
and the spline chain. 

According to the cubic approximation for S'(r), S'{t]^l) is a linear function of ^(rTVL-i)) 
S{tnl), S"{tnl-i) and S'^t^l), as shown by eq. (9). The "physical" equation (45) and 
the spline equation (9) lead to a new linear relation among S{tnl-i), S"{rNL-i), S{tnl) 
and S"{tnl) that, together with eq.s (39), (40) and (44) lead easily to the explicit values 
of the latter four variables. Consequently we easily obtain also the values of S'{tnl) and 
S"'{tnl)- The exphcit values of these variables at tnl allow us to compute those of the set 
{-'"■^(TiVL, A* j)} through eq. (5). 

Once the explicit values of S{tnl-i)-, S{tnl): S"{tnl-i): S"{tnl) as well as those of 
the set {I'^{tnl: pi'J: J — ^:^D)} are known, it is straightforward to compute those of the 
set {/+(TiVL-i), A«j, ^ = 1, ND)} via eq. (15). Then eq.s (27) and (28) will yield the explicit 
values of S{tnl-2) and S"{tnl-2), hence those of the set {I~^{tnl-2: I^J: J — ^D)}. And 
so on along the back- substitution. 

4. Conclusions 

Our Implicit Integral Method is based on the progressive treatment of the different 
layers that consitute a model of the stellar atmosphere physical system, from the outermost 
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layer (the surface) to the deepest one (the bottom) . The protagonist variables of the method 
are the upgoing and downgoing specific intensities /"""(r, fi) as well as the corresponding source 
functions that besides the thermal sources include a scattering-like integral into which there 
enter the foregoing specific intensities. Precisely, the study (and the elimination) of each 
single layer {tl,tl+i) leads to a relation that links linearly the value S{tl) of the source 
function at tl with S{tl+i), the value at tl+i- Once obtained via the study of the last layer 
the relation between the values of the source function at the two last optical depth points, 
the boundary condition at tatl given by eq. (5) makes it possible to compute the explicit 
values of S{tnl) and S{tnl-i), hence all the others. 

In order to design the required elimination scheme it is necessary to employ a math- 
ematical model for S{t). In principle the simplest and easiest model would be a piece- wise 
linear one, but the discontinuity of the first derivative S'{t) at each knot tl can imply severe 
errors and possible numerical instabilities because the above discontinuity is incompatible 
with the radiative transfer (RT) process itself, where both /"""(r, /x) and their first derivatives 
must be continuous, and therefore also the mean intensity J(t) and its first derivative. Thus 
the foregoing model cannot be correct, but for extreme cases of the thermal sources. 

A piece- wise parabolic model warrants the continuity of S'{t) at all depth points. 
Such a model shall include also S'(t) as a protagonist variable in the process of progressive 
elimination of the atmospheric layers. Hence 5"(r) must be put into relation with the 
foregoing protagonist variables, which can be done either mathematically or physically. 

Prom the mathematical standpoint we can introduce S'{t) by means of the formula 

S'{T,) = 2 ^(^^+0 - ^(^^) _ s'{r,^^) , (46) 
Tl+i — Tl 

which could however introduce numerical instabilities because of the difference between the 
two terms in the right-hand side, above all in the back-substitution process that works with 
explicit values. On the other hand, from the physical standpoint 5"(r) could be included 
by taking into account at all the optical depth points the equations (12) and (13) for the 
derivatives of the source function. However, in case that £(t) and B{t) show large variations 
(as it is the case of the formation of Lyman a in cool stars), severe instabilities may appear, 
too. 

These drawbacks can be avoided by introducing a piece-wise cubic approximation, 
where a further protagonis variable has to be included, namely the second derivative S"{t). 
That is, by means of a cubic spline model that automatically warrant the continuity of S{t) 
and its first two derivatives. To circumvent the explicit calculation of the derivatives makes 
the above difficulties vanish. 
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The source function S{tl) at each depth point tl will be expressed as a linear function 
of S{tl^i) and S"{tl^i). Therefore we shall transmit from any optical depth to the next one 
also the (implicit) value of S"{tl^i). This is achieved thanks to the fundamental relation 
that assures the continuity properties imposed by the cubic spline condition (cf. cq. [8]). 
However the propagation of S{t) and S"{t) via a cubic spline model constitutes a two- 
point boundary problem. Nevertheless this is perfectly compatible with the treatment of 
the transmission of the specific intensities, as it is performed in the scheme for the solution 
of the two-point boundary value RT problem. Both propagation processes can be treated 
simultaneously. 

Under these conditions wc can warrant the elimination of many of the causes of 
instability that can spoil the algorithm for the solution of the system of specific RT equations 
coupled through a scattering-like term in the source function, whose initial conditions are 
assigned at different points of the physical system. 
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